function [S, Rt, yss1] = KimShephardChib_Ivan3(ystar,h)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here

T = length(h);

% Mixture of Seven Normals (Table 4 from Kim, Shephard and Chib 1998)
% Note: means m_i already adjusted!!

q_i = [0.0073 .10556 .00002 .04395 .34001 .24566 .2575];
m_i = [-10.12999 -3.97281 -8.56686 2.77786 .61942 1.79518 -1.08819] - 1.2704;  %% means already adjusted!! %%
sig_i = [5.79596 2.61369 5.17950 .16735 .64009 .34023 1.26261];
sqrtsig_i = sqrt(sig_i);

Pr = repmat(q_i,T,1).*normpdf(repmat(ystar,1,7),2*repmat(h,1,7)+repmat(m_i,T,1), repmat(sqrtsig_i,T,1));
Pr = Pr./repmat(sum(Pr,2),1,7);

temprand = rand(T,1);
S = 7 - sum(repmat(temprand,1,7)<cumsum(Pr,2),2)+1;
S(S==8)=NaN; % 8 means that the observation is actually missing

% Construct Matrices for Kalman Filter
Rt = zeros(1,1,T);
yss1 = zeros(T,1);
for i = 1:T
    imix = S(i);
    if isfinite(imix)==1; Rt(1,1,i) = sig_i(imix); yss1(i,1) = ystar(i,1) - m_i(imix); 
    else  Rt(1,1,i) = NaN; yss1(i,1) = NaN; end
end

% [xtt, Sig, ~, ~,~] = KalmanFilter(yss1,[],[],h0,Sigh0,[],1,2,[],Q,Rt,[]);
% h = Bai_Wang2(1,Q,[],[],Sig,xtt,1);

end

